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What is the fuss all about? 
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The goal: 


Numerical Task = > Finite Computation ==> Distribution on Output 
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Optimisation 


x^ = arg max f(x) 


e Well-defined: 
e f € C%(&) for some a > 0 and x € X a compact subset of R7. 
ə Well-posed: 


ə Allowed n evaluations of f(-) at inputs which you can select. 
ə Aim to minimise ||X* — x* ||». 
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x* = arg max f(x) 


Two distinct requirements: 
e A method to select the function evaluation locations x1,..., Xn- 


e Uniform grid over A? 
ə Adaptive selection, e.g. gradient ascent with estimated gradients? 


e An estimator ((xi,f(xi)))i-i — X . 
ə The empirical maximum X€* = xj* where i* = arg maxj=1,....n f (xj)? 
e Something better? 


ə Key idea: Estimator uncertainty quantification! 
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Probabilistic Optimisation 


Start with the data that have been collected: 
D = {(xi, f(x)) Ha 


Bayesian linear regression onto a basis ($i): 
f(x) = &1éx(x) + +++ + Bmón(x) 


with n < m € NU {oo}. 
e Prior p(f1,.--; Bm) 
e Likelihood [T7 , 9(f(xi) — G1d1(xi) — --- — Bmómxi)) 
e Posterior p(fi,..., Bn|D) 


e Posterior marginal p(x*|D) 
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Probabilistic Optimisation 


Calculations for the conjugate set-up: 
CO (Lae (= (Gigs sesta) 
e Prior: BJA ~ N(O, A!) and A ~ p(A) « A. 


9 Posterior: 


B|D ~ MVT (u +8 p) Pf, a '$'(I«-99')"o9f)(I)- 99), n) 


where f = (f(x1),...,f(xn)) and [®]j = 9j(x:). 
ə Posterior marginal: x*|D ~ ? 
ə Draw 8 from B|D 


e Evaluate x* = arg max £1é1(x) +- -- + Ümóm(x) 
e Repeat. 
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Probabilistic Optimisation 


Compute x* — argmax f(x): 


Objective Function 
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Close connection between statistics and design of numerical optimisation methods. 


Similar to "Bayesian optimisation" (Mockus, 1989). 
Kernel trick maps to Gaussian processes. 
The distributional output p(x*|D) provides uncertainty quantification. 


Propagation and the Bayesian mantra of Dawid. 


Numerical analysts want to consider order of convergence and constants (of the 
point estimator). 


Similar considerations relevant to posterior contraction. 
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Application 


Objective function model 


Estimated objective function value 


https: //uk.mathworks.com/help/stats/ 
tune-random-forest-using-quantile-error-and-bayesian-optimization.html 
(but the MATLAB function doesn't provide uncertainty quantification!) 
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1= | fe 09s 


e Well-defined: 

e f € Lo(x)? 

e f € C? (X) for some a > 0 and x € X a compact subset of R. 
e Well-posed: 

ə Allowed n evaluations of f(-) at inputs which you can select. 

ə Aim to minimise |f — I|. 
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ə Key idea (again): Estimator uncertainty quantification! 


Chris. J. Oates Probabilistic Numerical Computation October 2017 @ Turing 


Probabilistic Integration 


Start with the data that have been collected: 


D = {(xi, f(x))Ha 


Chris. J. Oates Probabilistic Numerical Computation 


Probabilistic Integration 


Start with the data that have been collected: 
D = {(xi, f(x)) Ha 


Bayesian linear regression onto a basis ($i): 
f(x) = &1éx(x) + +++ + Bmón(x) 


with n < m € NU {oo}. 
e Prior p(£1,.... Bm) 
e Likelihood [[7 , ó(f(x;i) — 81¢1(x;) — +--+ — 8móm(xi)) 
e Posterior p(fi,..., Bn|D) 
e Posterior marginal p(/|D) 
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Probabilistic Integration 


Calculations for the conjugate set-up: 
O llas (= (Gig. se sea) 
ə Prior: BJA  N(0, M) and A ~ p(A) A1. 


9 Posterior: 


BID ~ MYT (u +8) o'f, -( '$'(I-$9')"9f)- 99), n) 


where f = (f(x1),...,f(xn)) and [E]; = ¢;(xi). 


ə Posterior marginal: 


jp = Student-T (W "(I+ & 9) *& f, 
EMO +06') '@f)v' (14+ 6'6) y, n) 


where [V]; = f ó;(x)n(x)dx. 
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Probabilistic Integration 


Compute f f(x)r(x)dx: 


Posterior Samples Corresponding Integrals 


OGD 


Reference Measure 
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ə Similar to "Bayesian quadrature’ (O'Hagan, 1991). 
@ Kernel trick maps to GPs. 
e Theoretical results were provided in Briol et a/. 2016: 
o Posterior mean Î satisfies / — | = Op(n-9/d*«). 
e Stronger assumptions on f, such as f € H%(0,1) ®--- & H%(0,1), lead to 
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o Select (xj)? , to minimise 


(x) w'(I4-995) 1v? 
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e Theoretical results were provided in Briol et a/. 2016: 


o Posterior mean Î satisfies / — | = Op(n-9/d*«). 
e Stronger assumptions on f, such as f € H%(0,1) ®--- & H%(0,1), lead to 
f — I = Op(n-9**) for an appropriately sparse basis {9i}. 
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e Posterior mean often coincides with standard quadrature methods (Särkkä et al., 
2016). 
@ How to select the (xi)—i? 
ə Recall the posterior scale was determined by 


1 
-(f'e'(I-ee')ef)v'(1-3- 6's)“ 
n 

o Select (xj)? , to minimise 


(x) w'(I4-995) 1? 


9 Select x; to minimise (*) based on Dir fixed? (Similar to "sequential Bayesian 
quadrature” .) 
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From Briol et al., 2016: 
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One of these is (a variant on) sequential Bayesian quadrature - which one? 


Application 


camera 


From Briol et a/., 2016. 
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Linear Algebra 


Ax =b 
e Well-defined: 
ə Aisa N x N symmetric positive definite matrix. 
e Well-posed: 
e Allowed n < N matrix-vector multiplications. 
ə Represented as [s; A|x = s! b for i = 1,..., n. 
e You are allowed to select the directions s;,..., Sp. 
ə Aim to minimise || — x||a where ||z||4 = Vz! Az. 
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e A method to select the directions 5s1,...,5;. 


e Random projections? 
e Sequential selection, e.g. gradient descent or conjugate gradient? 


e An estimator ((s; A, s; b)} 1 — &. 
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Linear Algebra 


Two distinct requirements: 
e A method to select the directions 5s1,...,5;. 


e Random projections? 
ə Sequential selection, e.g. gradient descent or conjugate gradient? 


o An estimator {(s; A, s; b))2., — &. 


ə A minimal ||- ||? norm vector that satisfies s! AX = s; b for i = 1,...,n? 
e Something better? 


9 Key idea (again): Estimator uncertainty quantification! 
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Recap: Conjugate Gradient Method 


Let (-,-)4 be the inner-product induced by || - |a. (i.e. (z,Z)a = z" AŽ.) 

Call z, € R" conjugate (w.r.t. A) if (z, Z)A = 0. 

Suppose that we have a conjugate basis {s1,..., sn} for R" (i.e. (sj, sj)A = O for all 
E) 


Consider the natural sequence of approximations 


n 
Xn = , QS; 
=L 


where each 


T (s.xX)a | sib 
j (Si, Si)A s As; 


can be computed in O(N?). The total computational cost is O(nV?). 
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Let (-,-)4 be the inner-product induced by || - |a. (i.e. (z,Z)a = z" AŽ.) 

Call z, € R" conjugate (w.r.t. A) if (z, Z)A = 0. 

Suppose that we have a conjugate basis {s1,..., sn} for R" (i.e. (sj, sj)A = O for all 
E) 


Consider the natural sequence of approximations 


n 
Ma = J QS; 
i=1 


where each 


T (s.xX)a | sib 
j (Si, Si)A s As; 


can be computed in O(N?). The total computational cost is O(nV?). 
So what is needed to proceed? 
e Need a smart choice of {51,..., sw]. 


e For theory, need to bound ||x — X,|| = || X oer a;x;|| in your favourite ||- ||. 
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Recap: Conjugate Gradient Method 


Aim is to adaptively select s; based on the computation up to iteration n — 1. 
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Aim is to adaptively select s; based on the computation up to iteration n — 1. 


Gradient Descent: Notice that x is a minimum of 


1 
(e) = 35x Ax — x b. 
This suggests to select s, = —Vf(Xn—1) which is equal to r;.1 = b — Ax; 1. However, 
this does not ensure {51,...,Sn} is a conjugate set. 


Chris. J. Oates Probabilistic Numerical Computation $ October 2017 @ Turing 


Recap: Conjugate Gradient Method 


Aim is to adaptively select s; based on the computation up to iteration n — 1. 


Gradient Descent: Notice that x is a minimum of 


HE) = Íx Ax - x b. 


This suggests to select Sn = —Vf(X, 1) which is equal to r; 1 = b — AX, 1. However, 
this does not ensure {51,...,5n} is a conjugate set. 


Conjugate Gradient: A more delicate procedure selects 


s; Arn—1 
Sn =  fn—i J Si 
r s.' As; 
icm i 


i.e. gradient descent plus Gram-Schmidt orthogonalisation w.r.t (-,-)a to subtract off 
components in the directions {s1,...,Sn—1} already used. 
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Aim is to adaptively select s; based on the computation up to iteration n — 1. 


Gradient Descent: Notice that x is a minimum of 


1 
HE) = 5x Ax — xb. 
This suggests to select Sn = —Vf(X, 1) which is equal to r; 1 = b — AX, 1. However, 
this does not ensure {51,...,5n} is a conjugate set. 


Conjugate Gradient: A more delicate procedure selects 


s; Arn—1 
Sn =  fn—i J Si 
r s.' As; 
icm i 


i.e. gradient descent plus Gram-Schmidt orthogonalisation w.r.t (-,-)a to subtract off 
components in the directions {s1,...,Sn—1} already used. 


For either method, the computational cost of selecting s; is O(N’), so the overall 
computational overhead added is O(nN?); the same order as random projections. 
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Probabilistic Linear Algebra 


Start with the data that have been collected: 


D = {(s; A, s; b))a 
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Probabilistic Linear Algebra 


Start with the data that have been collected: 


D = {(s; A, s; b))za 


Deploy full Bayesian inference for x: 
e Prior p(x) 
ə Likelihood [T^ , ó(s; Ax — s; b) 
9 Posterior p(x|D) 
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Probabilistic Linear Algebra 


Calculations for the conjugate set-up: 
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Calculations for the conjugate set-up: 
o bet. S= [irae] 


ə Prior: 
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Probabilistic Linear Algebra 


Calculations for the conjugate set-up: 
e Let: S — [sic esso] 
e Prior: x|\ ~ N(0, Al) and A ~ p(A) x AT?. 
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Probabilistic Linear Algebra 


Calculations for the conjugate set-up: 
e Let: S — [s.s] 
e Prior: x|\ ~ N(0, M) and À  p(A) x AW. 


ə Posterior: 
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Probabilistic Linear Algebra 


Calculations for the conjugate set-up: 
© Let: S = [sts] 
e Prior: x|\ ~ N(0, M) and À ~ p(A) x A7?. 


e Posterior: 


JD ~ MVT (( + A S SA)™*A" S" Sb, 


-(b S A' S + SAA' S) ^SASb)(I-- A' S' SA) ', n) 
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Probabilistic Linear Algebra 


Calculations for the conjugate set-up: 
© Let: S = [sts] 
e Prior: x|\ ~ N(0, M) and À ~ p(A) x A7?. 


e Posterior: 


xD MVT (( + A S SA) ‘A'S Sb, 


-(b S A' SQ + SAA' S) ^SASb)(I-- A' S' SA) ', n) 


ə This is for general S. 
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Probabilistic Linear Algebra 


Calculations for the conjugate set-up: 


O Lee S e busco" 
e Prior: x|\ ~ N(0, M) and À  p(A) x AW. 
9 Posterior: 


xD MVT (( + A S SA) ‘A'S Sb, 


-(b S A' SQ + SAA’ S) ^ SASb)(1 + A S SA) ^, n) 


ə This is for general S. 


ə For the conjugate gradient method applied to the pre-conditioned system 
A' Ax = A' b we have the orthogonality equation SAA'S' = I and the above can 
be further simplified. 
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Probabilistic Linear Algebra 


[no video for this one!] 
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Some Discussion Points 


e Approximate linear solvers used extensively in engineering applications. 
9 Also relevant in statistics, e.g. simulation of spatial random fields. 


e |t turns out that the posterior mean in our construction coincides with the classical 
conjugate gradient (CG) method applied to the pre-conditioned system 
A' Ax — A! b. 


ə Thus the classical error bounds for CG are inherited. 


e The full posterior can be computed in O(nN?). 
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From McAdams et al., SIGGRAPH 2010: 
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From McAdams et al., SIGGRAPH 2010: 
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Solution of Differential Equations 


E.g.: 


Au(x) 
u(x) 


F(x), XE 
0 xc0OX 


e Well-defined: 
o X C R7 be CH 
e f € LP(X), p» n/2 
ə Cor. 9.18 in Gilbarg and Trudinger ensures 3! solution u € W2^(x)n C?(x ux) 


loc 


Chris. J. Oates Probabilistic Numerical Computation 


Solution of Differential Equations 


E.g.: 
Au(x) = f(x), xexX 
u(x) = 0, xcoX 

e Well-defined: 

o X C R^ be CH 

o f € LP(X), p » n/2 

ə Cor. 9.18 in Gilbarg and Trudinger ensures 3! solution u € WeP(X) nc9(xuax) 
e Well-posed: 
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Solution of Differential Equations 


E.g.: 


Au(x) 
u(x) 


F(x), XE 
0 xc0OX 


e Well-defined: 
o X C R7 be CH 
e f € LP(X), p» n/2 
ə Cor. 9.18 in Gilbarg and Trudinger ensures 3! solution u € W2^(x) n C°(X U 8x) 


loc 


e Well-posed: 
e Allowed n evaluations of f(-) at inputs x1,..., x; € X U OX which you can select. 
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Solution of Differential Equations 


E.g.: 


Au(x) 
u(x) 


F(x), XE 
0 xc0OX 


e Well-defined: 
o X C R7 be CH 
e f € LP(X), p» n/2 
ə Cor. 9.18 in Gilbarg and Trudinger ensures 3! solution u € W2^(x)n C?(x ux) 


loc 


e Well-posed: 


e Allowed n evaluations of f(-) at inputs x1,..., x; € X U OX which you can select. 
ə Aim to minimise fy ||ô(x) — u(x)||2dx. 
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Au(x) = f(x), xexX 
u(x) = 0, xcoX 


Two distinct requirements: 
ə A method to select the function evaluation locations x1,...,Xn € X U OX. 


I = = = - 
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Two distinct requirements: 
ə A method to select the function evaluation locations x1,..., x». € X U OX. 


e Corners of a mesh on XY U OX? 
e An adaptive method? 
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Solution of Differential Equations 


Au(x) = f(x), xE 
u(x) = 0, xcoxXx 
Two distinct requirements: 
ə A method to select the function evaluation locations x1,...,Xn E€ X U OX. 


ə Corners of a mesh on ¥ U OX? 
e An adaptive method? 
e An estimator {(x;, f(x;i))] i-i  Q(-). 
ə Linear interpolation of the f(x;) and then solution of the PDE? 
e Something better? 


ə Key idea (again): Estimator uncertainty quantification! 
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Probabilistic Solution of Differential Equations 


Start with the data that have been collected: 


D = {(xi, f(x))Ha 
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Probabilistic Solution of Differential Equations 


Start with the data that have been collected: 
D = {(xi, f(x)) Ha 


Bayesian linear regression onto a basis ($i): 
f(x) = &1éx(x) + +++ + Bmón(x) 


with n € m € NU {oo}. 
e Prior p(£1,..., Bm) 
e Likelihood [T7 , 9(f(xi) — Gid1(xi) — --- — Bmóm(xi)) 
e Posterior p(fi,..., Bn|D) 
ə Posterior marginal p(u(-)|D) 
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Probabilistic Solution of Differential Equations 


Calculations for the conjugate set-up: 
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Probabilistic Solution of Differential Equations 


Calculations for the conjugate set-up: 


O [Es (= (Gipsa) 
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Probabilistic Solution of Differential Equations 


Calculations for the conjugate set-up: 


O [Ets (= (bin ees ole) 


9 Prior: 
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Calculations for the conjugate set-up: 


© llan (= (Gigs sse) 
e Prior: BJA ~ N(O, AJ) and A ~ p(A) A1. 
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Probabilistic Solution of Differential Equations 


Calculations for the conjugate set-up: 
© [bg (= (Gigs sse) 
e Prior: BJA ~ N(O, A!) and A ~ p(A) « A. 


9 Posterior: 
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where f = (f(x1),...,f(xn)) and [E]; = ¢;(xi). 
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Probabilistic Solution of Differential Equations 


Calculations for the conjugate set-up: 
© llan (= (Gigs sse) 
e Prior: BJA ~ N(O, M) and A ~ p(A) « AT. 
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where f = (f(»1),..., f(x;)) and [®] = ¢;(xi). 
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Probabilistic Solution of Differential Equations 


Calculations for the conjugate set-up: 
O [hae fe! (ens solem) 
e Prior: BJA  N(0, M) and A ~ p(A) A1. 


9 Posterior: 


BID ~ MYT (« +61) o'f, A '$'(I-$9')"9f)-9'9) , n) 


where f = (f(x1),...,f(xn)) and [9]; = ¢;(xi). 


ə Posterior marginal: 


u(x)D ~ Student-T (U(x) (1 +818) BTF, 
La t-99)"9f)U(x)' (1+ 9' 9) U(x), n) 


where [U(x)]; = u;(x) and u; solves Au = ¢; on X and u = 0 on OX. [N.B. Don't 
need to explicitly compute the à; if you have the Green's function of the PDE.] 
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Probabilistic Solution of Differential Equations 


Solve the ODE $ = f(x), u(0) = uo on x € [0, T]: 


Solution to an ODE 
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Some Discussion Points 


ə Posterior mean coincides with a classical “collocation” method. 


e Generalises to GPs with the kernel trick. 
e Theoretical results (for a method based on GPs) in Cockayne et al., 2016: 
e The posterior mean converges in ||- ||;5 at a rate 


O(h^-^- 2), 
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Some Discussion Points 


ə Posterior mean coincides with a classical “collocation” method. 


e Generalises to GPs with the kernel trick. 
e Theoretical results (for a method based on GPs) in Cockayne et al., 2016: 
e The posterior mean converges in || - ||;5 at a rate 


O(n*-?-$). 


ə The posterior mass for a ball of radius e centred on the true solution u(-) scales as 


20 —2p—d 
1-0 (——) 
€ 
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Application 


http://www.svflux.com/subdomains/svflux.com/index.shtml 
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Summary 


Theoretical Questions 


9 General theory? 


e Beyond linear and Gaussian assumptions? 


ə Experimental design? 


Chris. J. Oates 


Lots of work to do, but initial results in: 


Cockayne J, Oates CJ, Sullivan T, Girolami M 
Bayesian Probabilistic Numerical Methods 
arXiv:1702.03673 (2017) 


= = - = = - 
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Computational Questions 


e Propagation of uncertainty through a computational workflow? 


e Compatibility of multiple probabilistic numerical methods? 


14 sD ptapngt - 120182 +1 Opens e 
O qe Kame State App&caton = tna time [7287 "0 
4t How  dlistaomin ia Bow _coe_vedate LM 2011-09-18 102148 EX E 


[Fig: IBM High Performance Computation] 


The sophistication and scale of modern computer models creates an urgent need to better 
understand the propagation and accumulation of numerical error within arbitrary - often 
large - pipelines of computation, so that "numerical risk" to end-users can be controlled. 


Chris. J. Oates Probabilistic Numerical Computation October 2017 Q Turing 


